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Abstract 

We discuss a method of approximate model reduction for networks of bio¬ 
chemical reactions. This method can be applied to networks with polynomial 
or rational reaction rates and whose parameters are given by their orders of 
magnitude. In order to obtain reduced models we solve the problem of tropical 
equilibration that is a system of equations in max-plus algebra. In the case of 
networks with fast nonlinear cycles we have to compute the tropical equilibra¬ 
tions at least twice, once for the initial system and a second time for an extended 
system obtained by adding to the full system the differential equations satisfied 
by the conservation laws of the fast subsystem. Our method can be used for 
formal model reduction in computational systems biology. 
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1 Introduction 

Networks of chemical reactions are widely used in chemistry for modeling cataly¬ 
sis, combustion, chemical reactors, or in biology as models of signaling, metabolism, 
and gene regulation. In order to cope with growing amounts of data, these mod¬ 
els tend to be as comprehensive as possible by integrating many variables and 
processes with several different timescales. For most problems in computation 
and analysis of complex systems, the upper limit on the size of the system that 
can be studied has been reached. This limit can be very low, namely tens of 
variables for system identification, symbolic calculation or bifurcation of attrac¬ 
tors of large dynamical systems. Model reduction is a way to bypass these 
limitations by replacing large scale models with ones containing less parameters 
and variables, that are easier to analyse. 

There are several traditional numerical methods for reducing networks of 
chemical reactions. These methods, such as computational singular perturba¬ 
tion (CSP, [13), intrinsic low dimensional manifold (ILDM, ;T3]) exploit the 
separation of timescales of various processes and variables of the model. In 
dissipative systems, fast variables relax rapidly to some low dimensional attrac¬ 
tive manifold called invariant manifold [7J that carries the slow mode dynamics. 


1 


A projection of dynamical equations onto this manifold provides the reduced 
dynamics Ell¬ 
in the last decade, the rapidly growing field of computational and systems 
biology produced biochemical reactions networks models for cell physiology, 
of increasing size and complexity. Model reduction has been identified as a 
highest-priority challenge for these fields, expected to tame the complexity and 
simplify the analysis of biological models. However, these models came with 
peculiarities and the extant traditional model reduction methods are not en¬ 
tirely suitable for this endeavour. Biological models suffer from structural and 
parametric uncertainty and one rarely has precise information about kinetic pa¬ 
rameters. One of the main problem of computational biology is the parameter 
space exploration and analysis of possible model behaviors. Therefore, formal, 
symbolic, or semi-quantitative model reduction methods are more appropriate 
than numerical methods that need precise parameters. 

Formal model reduction can be based on conservation laws, exact lumping 
[B], and more generally, symmetry |3]. For chemical reactions networks with 
fast reaction cycles and fast species, lowest order approximations of attractive 
invariant manifolds are provided by quasi-equilibrium or quasi-steady state ap¬ 
proximations 0. These two approximations allow model reduction by graph 
reconstruction via linear lumping, pruning, and algebraic elimination of fast 
variables mm- Graphical reduction methods use elementary modes m , or 
the Laplacian defined on the graph of complexes of the reactions network m, 
but have little or no connection with singular perturbation methods and do 
not exploit multi-scaleness of biochemical networks. A fully formal reduction 
method exploiting orders of magnitude of variables and parameters is still miss¬ 
ing. 

In this paper we present a new model reduction method, inspired by tropical 
geometry and analysis. This method is particularly suited for computational 
biology because it combines graphical approaches, semi-quantitative reasoning 
and symbolic manipulation. 

The plan of the paper is the following. The second section introduces the 
biochemical reactions models and the tropical geometry concepts needed for 
the presentation of our results. In the third section we discuss the relation 
between tropical variety and Newton-Puiseux series. In the fourth section we 
provide general results of existence of an invariant manifold for biochemical 
systems with fast species and fast cycles. We also discuss how to choose slow 
and fast variables in this case, and how to define reduced models describing 
the slow dynamics on the invariant manifold and the fast relaxation towards 
this manifold. In the fifth section we apply these results to a nonlinear cycle of 
reactions. 


2 Definitions and settings 

We consider biochemical networks described by mass action kinetics 

( 2 . 1 ) 

j 

where kj > 0 are kinetic constants, Sij are the entries of the stoichiometric 
matrix (uniformly bounded integers, \Sij\ < s, s is small), cXj = (a{, a J 2 , ■ . ., a 3 n ) 
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are multi-indices, and x ai = x 1 1 x 2 2 .. .i„”. We consider that a\ are positive 
integers. However, the approach can be extended to rational or real indices. 

The choice is not restrictive, because most kinetic laws used in compu¬ 
tational biology can be decomposed into simpler steps each one obeying mass 
action law mug. Extensions of our approach, directly applicable to models 
whose rate functions are ratios of two polynomials (such as Michaelis-Menten 
or Hill functions) without expanding them into mass action elementary steps, 
were briefly discussed in m and will be presented in detail in future work. 
S-systems, used to model metabolic networks and for which ay are rational or 
real multi-indices |24j . are also covered by our approach. 

For slow/fast systems, the slow invariant manifold is approximated by a 
system of polynomial equations for the fast species. This crucial point allows us 
to find a connection with tropical geometry. We introduce now the terminology 
of tropical geometry needed for the presentation of our results, and refer to (15] 
for a comprehensive introduction to this field. 

Let /i, fi , • * •, fk be multivariate polynomials, fi £ C[xi,X 2 , ■ ■ ■ ,x n \. By 
considering sums of products of these polynomials by arbitrary polynomials we 
define the ideal I C C[xi, X 2 ,..., x n \ generated by them. The ideal is important 
in the context of solving systems of algebraic equations because any solution of 
the system fi(x) = 0, / 2 (a;) = 0,..., fk{x) = 0 is also a solution of /( x) = 0 
where f £ I. Important reasons for considering the generated ideal in the 
context of model reduction will be discussed in the Section [3] 

Let us now consider that variables Xi, i £ [l,n] are written as powers of 
a small positive parameter e, namely ay = aye ai , where ay has order zero. 
The orders cy indicate the order of magnitude of ay. Because e was chosen 
small, ay are lower for larger absolute values of ay. Furthermore, the order 
of magnitude of monomials x a is given by the dot product fj, =< a, a >, 
where a = (ay, 02 ,..., a n ). Again, smaller values of /x correspond to monomials 
with larger absolute values. For each multivariate polynomial / we define the 
tropical hypersurface T(f) as the set of vectors a £ R n such that the minimum 
of < a, a > over all monomials in / is attained for at least two monomials in 
/. In other words, / has at least two dominating monomials. 

A tropical prevariety is defined as the intersection of a finite number of 
tropical hypersurfaces, namely T(f 1 , / 2 , ■ ■■, fk) = ^ie[i,k]T(fi). 

A tropical variety is the intersection of all tropical hypersurfaces in the ideal 
/ generated by the polynomials / 1 , / 2 ,..., fk, namely T(I) = D/ 6 /T(/). The 
tropical variety is within the tropical prevariety, but the reciprocal property is 
not always true. 

For our purposes, we slightly modify the classical notion of tropical prevari¬ 
ety. A tropical equilibration is defined as a vector a £ R" such that < a, a > 
attains its minimum value at least twice for monomials of different signs, for 
each polynomial in the system / 1 , / 2) ■ • •, /fc- Thus, tropical equilibrations are 
subsets of the tropical prevariety. Our sign condition is needed because we 
are interested in approximating real positive solutions of polynomial systems 
(the sum of several dominant monomials of the same sign have no real strictly 
positive roots). 

In this paper we discuss how tropical equilibrations can be used for model 
reduction of chemical reactions networks. Tropical equilibrations indicate dom¬ 
inant monomials whose equality define approximated invariant manifolds. Fur- 
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thermore, they can be used to compute the timescales of the species, which is 
important for deciding which species are fast and can be eliminated by quasi- 
stationarity conditions. 

More precisely, we assume that parameters of the kinetic models (2.1) can 
be written as 

kj = kj£ J: ’. ( 2 - 2 ) 


The exponents 7 j are considered to be integer or rational. For instance, the 
following approximation produces integer exponents: 


7 .j =round(log(fcj)/log(e)), (2.3) 

where round stands for the closest integer (with half-integers rounded to even 
numbers). Without rounding to the closest integer, changing the parameter 
e should not introduce variations in the output of our method. Indeed, the 
tropical prevariety is independent on the choice of e. However, rounding to 
integer or rational exponents is needed in order to ensure that our lowest order 
approximations can be extended to series. 

Timescales of nonlinear systems depend not only on parameters but also on 
species concentrations, which are a priori unknown. We introduce the species 
orders vector a = (ai, 02 , • ■ •, a„), such that x = xe a . Of course, species orders 
vary in the concentration space and have to be calculated. To this aim, the 
network dynamics is first described by a rescaled ODE system 

^=Y^e' J :'- a 'k :l S tJ x^, (2.4) 

j 


where 


N = 7 3 + ( a > a 


■3 h 


(2.5) 


and (,) stands for the dot product. 


The r.h.s. of each equation in (2.4) is a sum of multivariate monomials in 


the concentrations. The orders /ij indicate how large are these monomials, in 
absolute value. A monomial of order pj dominates another monomial of order 
fly if Pj < Hj'- 


The timescale of a variable x.-, is given by T- ^ = A ^ whose order is: 


v% = min{/ij|5y ^ 0} - a. 


( 2 . 6 ) 


The order z/j indicates how fast is the variable 27 (if vp < Vj then xp is faster 
than v^) . 

The tropical equilibration problem consists in finding a species order vector 
a such that 

min (T j + (a, (*j)) = min hj + (“> a j)) ( 2 -7) 

0 J,Oij< 0 


We have shown in [IS] that species orders ai can be computed as solutions of 
(2.7). As discussed above, these solutions belong to the tropical prevariety of 
the polynomials defining the chemical kinetics. One of the problem of this ap¬ 
proach is that the tropical prevariety is too large, namely one can find too many 
tropical equilibrations. Although all these equilibrations can formally lead to 
reduced models, some are spurious. In this paper we propose to use the tropical 
equilibrations in a smaller set of the tropical variety. This choice is natural, 
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because by a result of Speyer and Sturmfels m the tropical variety is related 
to Newton-Puiseux series. As a matter of fact, choosing tropical solutions in the 
tropical variety ensures their lifting to series. Furthermore, using the tropical 
variety allows us to overcome another limitation of our previous application of 
tropical geometry ideas to model reduction. Namely, in |18| our reduced models 
were obtained by tropical truncation (consisting in neglecting dominated mono¬ 
mials and keeping only lowest order monomials in the differential equations). 
This method leads to unbounded errors when fast cycles are present in the re¬ 
actions network. Indeed, tropical truncation can generate fast subsystems that 
have conservation laws not present in the initial system. Although this kind 
of truncation is accurate on short timescales, it does not cope with slow relax¬ 
ation of the mass carried by the fast cycles. From a geometrical point of view, 
these conservation laws define sums of polynomials belonging to the ideal and 
contribute to the definition of the tropical variety. From a biochemical point of 
view, the conservation laws provide pools of species whose total mass relaxes 
slowly and should stand as supplementary slow variables. By this new approach, 
we use both pruning and pooling in order to reduce the biochemical reactions 
networks. 


3 Newton-Puiseux series and tropical equilibra¬ 
tions. 


In this section we introduce the Newton-Puiseux series and discuss their relation 
with the tropical equilibrations. 

By K = C((e 1 / 00 )) we denote the field of Newton-Puiseux series, i.e. all the 
series of the type 

x(e) = C\e~* + c 2 e^ + ..., (3.1) 

where Ci £ C, ai < a 2 < ... are integers, q is a positive integer. The series 
are convergent in some neighborhood of the origin, the origin being excluded if 
ai < 0. 

The Puiseux theorem [5] says that K is algebraically closed, i.e. that every 
nonconstant polynomial in K[x\ has a root in K. In particular, any root of 
a polynomial whose coefficients are powers of e can be written as a Newton- 
Puiseux series in e. Furthermore, the leading term order ^ can be calculated 
using the Newton polygon construction. Suppose we want to solve the equation 

P(x, = k j^ j x aj = 0, (3.2) 

i 


where 77 are integers and 07 are positive integers. In this case, Puiseux theorem 
ensures that Eq.(3.2) has solutions of the type (3.1). By substituting in (3.2) 
x(e) = C\e (1 + 27 (e)) (where 27 (e) denotes terms of order larger than zero in 
e) we get 

P{x, e) = q + ri(e) = 0, 

3 

where 77 (e) collects higher order terms. Necessary conditions for P(cc, e) = 0 
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read at lowest order 


E k i<? = 0 (3-3) 

j^j+aiocj /q=m 

m = min( 7 j + a±aj/q) (3.4) 


In order to satisfy (3.3), the minimum in (3.4) should be attained at least twice. 


Furthermore, if one looks for real solutions c\ £ 
that at least t\ 
signs, namely: 


then from (3.3) it follows 
that at least two kj corresponding to the minimum (3.4) should have opposite 


minj >kj> o('Yj + ai atj/q) = + ah atj/q). 


(3.5) 


We should note that (3.5) is a necessary, but not sufficient condition for real solu¬ 
tions (for instance x 2 —x+l = 0 satisfies the condition but has no real solutions). 
The above condition means that the lowest order a\/q in the Newton-Puiseux 
series solution has to satisfy a tropical equilibration problem. Geometrically, 
— a\/q is the slope of an edge of the Newton polygon, defined as the upper 
convex hull of the points of planar coordinates (ay, 7 j) (i.e. the convex hull 
including with any point the vertical half-line emanating up from this point). 
For instance, the leading terms in solutions of x 3 + ex 2 — x + e 2 =0 have orders 
e° or e 2 (see Figure [l]). The Newton polygon method can be generalized to mul¬ 
tivariate polynomials and we have implemented it in an automatic algorithm 
for computing tropical equilibrations presented elsewhere |2.3| . 



a) 


Figure 1: Newton polygon and roots of the polynomial x 3 + ex 2 — x + e 2 . a) The 
Newton polygon edges are indicated by thick lines and the limiting monomials all 
satisfy the sign condition. The slopes of the edges are 0 and —2 corresponding to 
leading terms in the Newton-Puiseux series of orders e° or e 2 , respectively, b) The 
absolute values of the polynomial’s roots are represented vs. e, in logarithmic scale; 
the slopes are the roots valuations. 

Fast variables of chemical reactions networks with multiple timescales satisfy 
quasi-stationarity equations. These are multivariate polynomial equations of the 
type 

P(xi,x 2 ,...,x n ,e) = E fc j e7i 4 ai)l 4 a ' )2 • ■ ■x i n j)ri = 0, (3.6) 

3 
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where 7 j are integers and ( aj)k are positive integers. Let us note that P(x) £ 
K[x\,X 2 , ■ ■ ■ ,x n \. Like in the case of the univariate equation (3.2), the tropical 


equilibrations provide lowest order approximations of the solutions of (3.6). We 


want to represent solutions of (3.6) by series, in other words we want to lift the 


tropical solutions to Newton-Puiseux series. In the univariate case, this is always 
possible by the Puiseux theorem. In the multivariate case, the lifting as Newton- 
Puiseux series of any tropical solution is ensured by a theorem of Kapranov [3] . 
In order to formulate this result let us introduce the valuation function Val(x) 
defined as the lowest power of e occurring in some Newton-Puiseux series x(e). 
When applied to species concentrations Xi , kinetic parameters A 7 and monomials 
kiX ai the valuation gives the orders a^, 7 j and /ij, respectively. As it can be 
easily checked 

Val(x(e)) = limlog e |a;(e)| 

€—>•0 

Suppose that the tropical equilibration condition defined in the introduction 


(Eq.(2.7)) is satisfied. By the same method as in the univariate case it can be 


shown that (2.7) is a necessary condition to have real Newton-Puiseux solutions. 


The Kapranov theorem 53 states that (2.7) is also a sufficient condition for 
having Newton-Puiseux solutions with prescribed lowest orders (ai, 0 , 2 ,.. •, a n ). 
More precisely, there are 27 (e) £ K such that V al(xi) = at and such that 
P(x 1 (e),x 2 (e),.. ,,x n (e),e) = 0. 

There is no analogue of Kapranov theorem working for systems of equa¬ 
tions. In this case, the condition (2.7) though necessary, is not sufficient for 


guaranteeing the existence of Newton-Puiseux solutions. 

In general, in order to obtain tropical equilibrations that can be lifted to 
Newton-Puiseux series we need to find a so-called tropical basis [lj. Let us 
consider that we want to find approximate solutions of n equations of the form 


(3.6), namely 


Pi{x, e) = 0, P 2 {x,e) = 0, ... ,P n (x,e) = 0 


(3.7) 


We first look for vectors (ai, a 2 , ■ ■ ■, a„) € K" satisfying the tropical equilibra¬ 
tion condition (2.7) for each polynomial Pfc, k £ [1 ,n]. This set is included in 
the tropical prevariety. Contrary to the case of one equation, it is not longer 
guaranteed that there are solutions Xi(e),x 2 (e),... ,x n (e) of (3.7) such that 


K dl{x 1, X 2 , . . . , 2? n ) — (on, C&23 • ■ * ) ^n)i 


(3.8) 


where Val means here application of the valuation componentwise. Solutions 


of (3.7) that satisfy (3.8) can be found if we solve a more complex problem. 


Let us supplement the system (3.7) with sums of products of the polynomi¬ 
als Pi, P 2 ,..., P n by arbitrary polynomials and solve the tropical equilibration 
problem for the augmented system. In other words, let us look for solutions in 
the tropical variety of the ideal / generated by P\,P 2 ,... ,P n - By a result of m 
in this case there are Newton-Puiseux solutions that satisfy the property (3.8). 


Although the ideal has an infinite number of elements, it can be shown that 
it is enough to solve the tropical problem for a finite set of polynomials in the 
ideal. A finite set of polynomials fi, f 2 , ■ ■ ■, ft generating the ideal / and such 
that T(/i) nr(/ 2 )n...nr(/ t ) = T(I) is called tropical basis. An algorithm for 
computing a tropical basis can be found in |Tj. However, the complexity of this 
algorithm can be double-exponential in the size of the system, both in time and 
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in space. In the remaining of this section we state a simple method for finding 
tropical solutions that can be lifted to Newton-Puiseux series. 

Generically, a system of n tropical equations in n variables has a finite num¬ 
ber of solutions. Indeed, the equality of two n-variate monomials corresponds to 
a (n — l)-dimensional hyperplane in the space of coordinates log (a;*), or, equiv¬ 
alently, in the orders a^. The intersection of n hyperplanes of dimension (n — 1) 
in a n-dimensional space is generically a point. Because the combinations of 
monomials that can be equilibrated are finite, the total number of solutions 
is finite. However, chemical reactions systems often have infinite branches of 
tropical equilibrations. For instance, a reversible reaction can equilibrate both 
reactants and reaction products. In quasi-equilibrium conditions, the forward 
and reverse rate monomials of this reaction are dominant, have equal orders 
and occur in equilibration equations of several variables. Therefore, some of the 
hyperplanes resulting from different tropical equations coincide. In these cases, 
there are infinite sets of tropical equilibrations. 

As an example, we can consider the following system of equations: 


y — x — ex 4 = 0 
x — y + ey 2 = 0 


(3.9) 


The valuations ai = Val(x(e)) and a 2 = Val(y(e)) of solutions of (3.91 satisfy 
the tropical equations min(ai, 1 + 4ai) = min(ai, 1 + 2a 2 ) = a 2 . The condition 
min(ai,l + 4ai) = ai, leads to the infinite branch of solutions ai = a 2 > 
— 1/3. The condition min(ai, 1 + 4ai) = 1 + 4ai leads to min(oi, 1 + 2a 2 ) = 
min(ai,3 + 8ai) = 3 + 8ai = 1 + 4oi, therefore a\ = —1/2, a 2 = —1, or 
to min(ai,3 + 8ai) = ai = 1 + 4ai, hence ai = a 2 = —1/3 solution already 


found. Thus the system (3.91 has an infinite branch of tropical equilbrations 


ai = a 2 > —1/3 and a isolated solution ai = —1/2, a 2 = —1. 

However, only two of these solutions lead to Newton-Puiseux solutions of 


(3.9). Indeed, the system (3.9) has 7 complex solutions, namely (0,0), (x, ±x 2 


where x is a solution of ex'^fx +1 = 0. Using the Newton polygon construction 
we find that the possible valuations of x are 0 or —1/2. It follows that valuations 
of real Newton-Puiseux solutions of Eq. (3.9) are (0,0), or (—1/2,—1). The 
only tropical solutions leading to Newton-Puiseux series is (0, 0), a point on the 
continuous, infinite branch of tropical solutions and (—1/2, —1), the isolated 
solution. 

We conjecture that all the isolated tropical equilibrations can be extended 
to Newton-Puiseux series. Therefore, if by supplementing the original system 
with sums or products of the original equations with arbitrary polynomials (i.e. 
considering the ideal generated by these equations) we find a system with only 
isolated tropical equilibrations, we believe that all of them can be lifted to 
Newton-Puiseux series. For instance, in the above example, let us consider to¬ 


gether to the equations (3.9) also their sum e(y 2 —x 4 ) = 0 and solve the resulting 
extended tropical system min(ai, 1 +4ai) = min(ai, 1 + 2a 2 ) = a 2 , ai = a 2 /2. 
This system has only two solutions, (0, 0) and (—1/2, —1). These two solutions 
are isolated. We have already shown that they can be lifted to Newton-Puiseux 
series. An ansatz for finding linear combinations of equations leading to isolated 
tropical equilibrations is to consider conservations laws of the fast subsystem. 
This ansatz will be used in the Sections 4|5 









4 Model reduction of biochemical reaction newt- 
works with fast cycles. 


In this section we introduce our model reduction method. We also state and 
prove our main result on the existence of invariant manifolds for networks with 
fast species and fast cycles. 

Let us call tropically truncated system associated to the tropical equilibration 
(ai, a 2 ,,a n ), the system obtained by pruning all the dominated monomials 
of ( |2.1| revealed by the rescaling ( |2.4[ ), i.e. 

^ = H *e[l,n] (4.1) 

teJ(i) 


where J(i ) = argmin^-, Sjj y 0) is the set of dominating reaction rates of 


reactions acting on species i and Hj are defined by (2.51. 

Like in the introduction, we introduce the orders = pj^ — ai, with pj^ = 
min (pj, 0). The rescaled truncated system reads 


dxj 

d t 


e Ui kjSijX aj , i£ [l,n]. 


(4.2) 


Variables Xi with smaller orders y t are faster. After reordering the variables we 
can consider that V\ < z/ 2 < ... < v n . Let us assume that the following gap 
condition is fulfilled: 


there is / < n such that Vf+i — Vf > 0, 


(4.3) 


meaning that two groups of variables have separated timescales. The variables 
Xf = (xi, X 2 , ■ ■ ■, Xf) are fast (change significantly on timescales of order of 
magnitude or shorter). The remaining variables x s = (xf + i, Xf + %, ..., x n ) 
are slow (have little variation on timescales of order of magnitude e~ u f). 

We call linear conservation law of a system of differential equations, a linear 
form c{x) =< c, x >= ciaq + C2X2 + ... + c n x n that is identically constant on 
trajectories of the system. It can be easily checked that vectors in the left kernel 
Ker l (S) of the stoichiometric matrix S provide linear conservation laws of the 


system (2.1l. Indeed, system (2.1) reads ^ = SR(x), where the components 


of the vector R are Rj{x) = kjX aj . If cr S = 0, then 
where c T = (ci, c 2 , ..., c n ). 


d<C,£C> 

d£ 


= c 1 SR(x) = 0, 


Let us assume that the truncated system (4.1), restricted to the fast vari¬ 
ables has a number of independent, linear conservation laws, defined by the 
vectors C\, c 2 , ..., c r , where Cfc = (cfci, Cfc 2 , ..., Cfc ra ). These conservation laws 
can be calculated by recasting the truncated system as the product of a new 
stoichiometric matrix and a vector of monomial rate functions and further com¬ 
puting left kernel vectors of the new stoichiometric matrix. We further assume 


that the truncated system and the full system (4.1) have no conservation laws 
in common. 

We now define the new variables yk = c kiXi, where k £ [1, r). These new 
variables satisfy the equations 




(4.4) 
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Let (oq, 02, ■ • •, a n ) be a solution of the tropical equilibration problem for the 
augmented system obtained by putting together (2.1) and (4.4). We define b k = 


min (ai,c k i 0) and p k = y Jc (k) ~ bk where p-j c ( k ) = min (p,j,c k i ± 0, Sy ^ 0), 

J c {k) = argrnin (pj,c k i 0, Sy 7^ 0). In rescaled variables y k = y k e bk we have 

3 

the following truncated rescaled system 


d y k 
d t 


z pk E E c ^^- 

jeJ c {k) 1 


(4.5) 


We assume that Vf < p k , meaning that the variables y k ,k £ [l,r] are slower 
than the variables Xi,i £ [1,7']. 

Since we have r conservation laws, we can eliminate r fast variables from the 
truncated system. One can suppose that these fast variables are £/_ r +1, ...,£/ 
can be expressed via the remaining variables a 'i, i £ [1, / — r) and y k , k £ [1, r]. 
Let us introduce vectors 


X r — (xi, X2, ■ • •, xj_ r ), X s — ('£f+ 1, •£/+2, • • •, x n ), y — (yi, 2/2, ■ • •, 2/r)* 

Then any function of x can be expressed via X ri X s and y. As a result, we 
obtain the following decomposition 


d Xi 

d t 


= e v *F i {X r ,y,X s ,e) i e [1 ,/-r], 


(4.6) 


where 


F i (X r ,y,X s ,e) = ^ k 3 S t 
a e J(i) 


x a *, 


daq 
d t 


= e I '>S i (X r ,y,X s ,e), i£[f+l,n], 


Si(X r ,y,X s ,e)= 

iG J{i) 

d y k 


d t 


= t Pk Y k (X r , y, X s , e), k£[l,r], 


(4.7) 


(4.8) 


Y k (X r ,y 1 X s ,e)= ^ £c w Sy®°y, 

jeJc(fc) 1 

where Y. S and F are analytic functions. 

The system (4.6) describes the evolution of fast modes. Because it was 
obtained from the truncated versions of the first f — r equations of (2.1), let 
us call it the fast truncated subsystem. As a matter of fact, the system (4.6) 
coincides with the first / — r equations of (4.2). 

Let us recall some notions of the dynamical systems theory. Let dx/dt = 
F{x) be a system of ordinary differential equations defined on an open domain 
FI of an Euclidean space with smooth boundary. Here F is a smooth function, 
for example, F £ C r , where r > 1. Let us consider an equilibrium (steady 
state) 0 (i.e., the relation F(cf>) = 0 holds) of this system. Let A be a linear 
operator that gives a linearization of r.h.s. of this system at 0: 


F(x) = A(x — 0) + 0{{x — 0) 2 ). 


10 














We say that this equilibrium is hyperbolic pS], if the distance d between the 
spectrum Spec a of A and the imaginary axis / = {z £ C : Re z = 0} is not 
zero: 


d = dist{SpecA , I) f 0. 


(4.9) 


If the spectrum lies in the left-half plane {z £ C : Re z < 0}, then this equi¬ 
librium is stable and locally attracting. In our case all systems depend on the 
parameter e > 0, therefore, d in (4.9) can depend on e. 

We can now formulate our main result. 


Theorem 4.1 Assume the gap condition (4.3) holds and that Vf < p^, k £ 
[1,7-]. Assume that for all values y and X s the fast truncated system (4.6) has 
a stable hyperbolic steady state 

Xi = <j>i(y, X a ) + higher order terms , i £ [1, / — r], 

such that the distance d(e) admits the estimate 

d > C 0 e K 

where k > 0 is small enough and Co is independent on e. Then for sufficiently 
small e > 0 system (4.1) has a locally attracting and locally invariant normally 
hyperbolic C p (p > 1) smooth manifold defined by 

Xi = (f>i(y, X s , e) + higher order terms , i £ [1, / — r], (4-10) 

and the dynamics of the slow variables y,Xf + i,Xf + 2 ,... ,x n for large times takes 
the form 

d Xi 


—= e Vi kjSijX + higher order terms , i £ [f + l,n] 

jeJ(f+ 1 ) 


df 

d fjk 
df 


cPk 


E E cuSijxf 3 + higher order terms, k £ [l,r] (4.11) 


jGJc(k) i 


where xf 3 = 


^2 2 


J f ^f-\-l x /+2 • • • 


Remark: “higher order terms” means some smooth functions of y,X s ,e 
having higher orders in e with respect to principal terms. 

Proof. We use the standard result, which follows from the theory of normally 
hyperbolic invariant manifolds m Consider the system of differential equations 


du/dr = Au + A F(u, v , A) + H(u, v, A), 
dv/dr = X^S(u, v, A), p > 0, 


(4.12) 

(4.13) 


where u £ K", v £ R m , A is a linear operator such that the spectrum of A 
satisfies (4.9) and lies in the left-half plane, F,G and H are smooth functions 
uniformly bounded in C k -norm on x R m x [0,1] for some k > 1. Moreover, 
H = 0(|w| 2 ) as u —y 0. 

It is clear that this system becomes slow/fast for small A > 0, where u 
are fast and v are slow. For any p < k and sufficiently small A there exists 
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a normally hyperbolic smooth locally attracting invariant manifold close to 0 : 
u = XU(v, A), where U is bounded in C p -norm. 


To apply this result, we reduce our system (2.1l to the form (4.12), (4.13). 


We introduce u = X r — 4> and makes a time change r = e K t. We introduce 
the variables v by v = ( X s ,y). Then, if k > 0 is small enough we obtain that 


system (4.6), (4.7), (|4.8|) can be rewritten in the form (4.12), (4.13) with A = e p 


for some p > 0. This completes the proof. 


5 A simple nonlinear cycle example. 

Let us consider the following example of a cycle of reactions that includes a 
complex formation reaction: 


a 1 ^a 2 ^a 3 ^a 1 , 


All 


kd 

a 2 ^ a 3 

fcs 


The mass action chemical kinetic equations for this cycle read: 

-k\X\ + k 3 X 3 — k 4 X 3 X 2 + k 5 x 3 

k 3 xi — k 2 x 2 - kiX\x 2 + k 3 x 3 
— k 3 X 3 + k 2 X 2 + k±X\X 2 — k 5 x 3 


daq 
d t 
dx 2 
d t 
dx 3 
d t 


(5.1) 

(5.2) 

(5.3) 


Consider kinetic constants that scale like ki ~ e Ti . For instance if e = 1/10 and 
k\ =l,k 2 = 0 . 1 , k 3 = 0 . 01 , = 0 . 01 , k 3 = 0.001 we get 


7i = 0 ,72 = 1 ,73 = 74 = 2 ,75 = 3. 


(5.4) 


The tropical equilibrations for this model are the solutions of the following 
min-plus equations 

min(7i+ai, 74 + 01 + 02 ) = min( 7 3 , 7 5 )+a 3 = min(7 2 +a 2 , 74 + 01 + 02 ) = min( 7 i+ai, 7 5 +a 3 ), 

(5.5) 


where the first, second and third equality of (5.5) follow from (5.1), (5.2), (5.3), 


respectively. Because q 3 < 75 we have min( 713 , 75 ) = q 3 . From (5.5) it follows 


that q 3 +a 3 = min (71 + 04,75 + a 3 ). Furthermore min( 7 i +ai, 75 +a 3 ) = 71+01 
(because 75+03 > 7 3 +a 3 ). Hence, in this case, the system of min-plus equations 
can be simplified to 

min( 7 i + oi ,74 + ai+a 2 ) = 7 3 + a 3 = min( 7 2 + a 2 , 7 4 + ai+a 2 ) = 71+04 (5.6) 


Only one of the possible outputs of the first min operation (5.6) has to be 


considered, namely min( 7 i + Oi, 74 + 01 + a 2 ) = 71 + ai, whereas the second min 
leads to two situations. It follows that there are thus two branches of tropical 
solutions, namely: 


ai> 7 2 - 74 , a 2 =01+71-72, 03 = 01 + 71-73 

ai < 72 - 74, a 2 = 71 - 74, «3 = ai + 71 - 73 


(5.7) 

(5.8) 


The tropical equilibration solutions vary continuously on each of the branch and 
are non-isolated. It is thus possible that some tropical solutions or an entire 
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branch can not be lifted to Newton-Puiseux series. In order to find solutions 


(5.31, 


that can be lifted we will consider equations in the ideal of (5.1 1 , (5.2 1 
starting with the conservation laws of the fast subsystem. It will come out 
that this is enough, as adding these conservation laws leads to isolated tropical 
equilibrations. 


Let us first consider the branch (5.71. By keeping dominating monomials 


of lowest order in e and pruning all the others we get the following truncated 
system: 

e 71 (—£727 + k 3 x 3 ) 

e' y2 (k 1 Xi - k 2 x 2 ) (5.9) 

e 73 {—k 3 x 3 + k 2 x 2 ), 


(T) 


dtci 
d t 

= e 71 

ax 2 
d t 

= e 72 

dx 3 

— f 73 

d t 



The truncated system (5.91 have the linear first integral (conservation law) 

y = xi+x 2 + x 3 . (5.10) 

The variable y is not a first integral of the full system, which implies that the 
truncated system (T) can not be a good approximation at large times. The 
exact dynamics of y is obtained by summing the equations (5.1 1 ,(5.21,(5.3): 


d V , 

— = -kiXix 2 

at 


k 5 x 3 . 


( 5 .H) 


We consider that y ~ e a “ and further equilibrate the equations (5.10),( |5TT l. 
We therefore get two more min-plus equations: 




* 2 / = min(ai,a 2 ,a 3 ) 
74 + aq + a 2 = 75 + a 3 


(5.12) 

(5.13) 


Assume the particular choice (5.4) of parameter orders. Then, for the trop¬ 
ical solution (5.7), it follows ai = 0, a 2 = —1, a y = a 3 = —2, zq = 0, v 2 = 1, 
v 3 = 2, u y = 3 which means that y is slower than Xi, 1 < i < 3. The resulting 
tropical equilibration is in fact unique and thus isolated. Indeed, considering 


the second branch of tropical equilibrations for the variables 27 , 27 , 2:3 (5.8) we 


find that y can not be equilibrated because (5.131 and (5.8) imply 75 = 73 which 


is not satisfied. The polynomial defining the dynamics of y being in the ideal of 


polynomials defining the dynamics of 27 , 3 : 2 , 2 : 3 , it follows that the branch (5.8) 
is not in the tropical variety and can be discarded. 

We will now use this isolated tropical equilibration to obtain reduced models. 


We will discuss three reduced models: the truncated model (T) in (5.9) that 
describes the relaxation dynamics towards the attractive invariant manifold, 


the reduced model (4.11) given by Theorem 4.1 and describing dynamics on the 


invariant manifold, and a third reduced model combining these two. 


The truncated system (5.9) copes only with the fast relaxation onto the 
invariant manifold. The tropical approximation of the invariant manifold is 


obtained by setting the l.h.s of (5.9) to zero, i.e. computing the steady states of 
the truncated system. This approximation is the half-line x 2 = k 3 k 2 1 xi, x 3 = 
kik^ 1 Xi, 27 < k 2 k^ x . By using the new tropically truncated equation: 


V = x 3 , 


we compute 27 , 27 , 27 from y and obtain the reduced model: 


dy ^ —j Ej. 3 Z- -1_ Ire'll 


J.2 


(5.14) 

























The reduced model ( R ) (5.15) copes with the slow dynamics on the invariant 
manifold. 

In this particular example the two approximations (T) and ( R ) are compos- 
able, i.e. they can be merged in a model with broader validity. By replacing 


y with x 3 in (5.15) and combining the resulting equations with the truncated 


system (T) we get the following model: 



(M) % = 


-hx 1 + k 3 x 3 
kixi — fc 2 x 2 
—k 3 x 3 + fc 2 x 2 — 


l k 2 1 k 3 kiX 3 + k 3 x 3 


(5.16) 


The model (M) is a multiscale reduction as it gives accurate approximations 
of both fast and slow parts of the trajectories. 

The comparison of different approximations is shown in Figure [2] The va¬ 
lidity of the multiscale reduction could depend on the initial data. We have 
determined numerically the domain of initial data leading to accurate reduc¬ 
tions. Starting from the same initial data we have integrated the full model 
(5.1),(5.2),(5.3) and the reduced model ( |5.16 l and obtained the trajectories x(t) 
and xr(t ), respectively. We have computed the error such as Hausdorff-Pompeiu 
distance between the sets {( logio(t ), logio(xi(t)),logio(x 2 (t)), logio(x 3 (t))), variable t} 
and {(logio(t),logio(xri(t)),logio(xr 2 (t)),logio(xr 3 (t))),variab\et}. We notice 
in Figure [3^t) that we can change the initial data on 7 decades and still keep 
the trajectories of the reduced model ( 5.16| ) close to the trajectories of the full 
model (|5T|),(p|),(pl). 


D 



Figure 2: Two trajectories of the nonlinear cycle model defined by 
Eqs.(5.1),(5.2),(5.3) starting from 0\ and 0 2 are represented in red solid line. 
The blue circles are the trajectories starting from 0\ and 0 2 computed with 
the tropical truncated model (T) defined by Eqs.(5.9). The blue crosses are 
the trajectories computed with the reduced model defined by Eqs.(5.15). A is 
the stable steady state of the model. The half-lines BC and BD belong to 
parts of tropical variety corresponding to the tropical solutions a\ > — 1 , a 2 = 
01 — 1 , a .3 = ai — 2 and a\ < — 1 , a 2 = — 2 , 03 = a\ — 2 , respectively. 
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Figure 3: Trajectories of full (Eqs.(5.11,(5.2),(5.3>) and multiscale reduced 
model (Eqs.(5.16)) were computed starting from the same initial condi¬ 
tions and the Hausdorff-Pompeiu distance between the corresponding sets 
{(log w (t),log w {xi), log 10 (x 2 ), log 10 (x 3 )) 1 variable t} were calculated. The ini¬ 
tial conditions where chosen from an uniform grid in logarithmic scale —5 < 
logio(xi) < 4, —5 < logio(x 2 ) < 2, —5 < logio{xz) < 4. a) The positions of 
the initial data leading to Hausdorff-Pompeiu distance less than 0.15 are shown 
by circles with color coded values of this distance. The lines indicates the same 
parts of tropical variety as in Fig[2] Initial data can vary on 7 decades with 
global relative errors less than 1 — 10 015 ~ 0.4 which for e = 1/10 and stiff 
trajectories is remarkably robust, b) The distribution of Hausdorff-Pompeiu 
distances is shown for the set of initial data. 


6 Conclusion 

We have shown how to relate the tropical equilibration problem to the slow- 
fast decomposition and model reduction of biochemical reactions networks. In 
the case of biochemical networks with mass action kinetics, we use tropical 
equilibration solutions to find which species are fast and which are slow. We have 
proposed elsewhere two methods for solving the tropical equilibration problem, 
a first one by reformulating it as a constraint satisfaction problem [26) and a 
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second one based on the Newton polyhedron [231 . 

Under rather general conditions, existence of small dimensional attractive 
invariant manifolds for reactions networks with fast cycles and species is shown. 

Our model reduction recipe consists in calculating tropical equilibration solu¬ 
tions at least twice. At the first step we solve the tropical equilibration problem 
for the initial system of differential equations. This allows us to identify the fast 
species, that constitute the fast subsystem of the model. The fast truncated sys¬ 
tem, obtained by pruning dominated monomials in the equations of fast species, 
copes with relaxation towards the attracting invariant manifold. The tropical 
equilibrations calculated at this step belong to the tropical prevariety, but not 
all of them lead to a reduced model. In order to filter tropical solutions we 
solve the tropical equilibration problem a second time. If the fast truncated 
subsystem has conservation laws different from the ones of the full system, we 
use them to define new slow variables. At the second step, we solve the tropical 
equilibration problem for the augmented system that is obtained by adding to 
the initial system the differential equations satisfied by the conservation laws 
of the fast truncated subsystem. These new equations are linear combinations 
of the initial ones. We conjecture that if the resulting solutions are isolated, 
then they belong to the tropical variety. Also, they lead to reduced models ob¬ 
tained by expressing fast variables in terms of the slow variables. The resulting 
reduced model copes with the dynamics on the invariant manifold. If after the 
second step one still has truncated systems with conservation laws and continu¬ 
ous branches of tropical equilibrations, the first two steps can be reiterated until 
there are no conservation laws different from the ones of the full model and all 
tropical equilibrations are isolated. 

Our method can be used as recipe for formal model reduction in compu¬ 
tational biology. Some steps of the recipe are already automated, such as the 
calculation of tropical equilibrations. The computation of conservation laws 
can be performed by methods from |25j . but we don’t exclude the existence of 
difficult cases when conservation laws are not enough for grasping the tropical 
variety. In these difficult cases, calculation of tropical basis can use methods 
from pQ. Another challenging step is the elimination of fast variables as solu¬ 
tions of systems of polynomial equations. When the polynomials of the fast 
truncated system contain only two monomials, we can apply rapid methods for 
toric systems mm- In general, fast truncated systems are fewnomials (have a 
small number of monomials) and can be approached by the methods for sparse 
polynomial systems m- However, even with fewnomials, there could be models 
for which the elimination of fast variables is difficult. Numerical methods can 
be used, as a last resort, to solve problematic cases. 

For a simple example, we suggested, without providing a general recipe, how 
to combine one-scale approximations to obtain a multiscale approximation that 
is valid on both fast and slow time scales. A general method for obtaining 
multiscale approximations was given in [5j for networks of monomolecular reac¬ 
tions (in these networks each reaction has at most one reactant and at most one 
product and the reaction rates are given by the mass action law) with separated 
kinetic constants. Multiscale approximations of nonlinear networks are much 
more difficult and will be discussed elsewhere. 
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